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[t] Glacierized change in the Himalayas affects river- 
discharge, hydro-energy and agricultural production, and 
Glacial Lake Outburst Flood potential, but its quantification 
and extent of impacts remains highly uncertain. Here we 
present conservative, comprehensive and quantitative pre- 
dictions for glacier area and meltwater flux changes in 
Bhutan, monsoonal Himalayas. In particular, we quantify 
the uncertainties associated with the glacier area and melt- 
water flux changes due to uncertainty in climate data, a 
critical problem for much of High Asia. Based on a suite 
of gridded climate data and a robust glacier melt model, 
our results show that glacier area and meltwater change 
projections can vary by an order of magnitude for different 
climate datasets. However, the most conservative results 
indicate that, even if climate were to remain at the present- 
day mean values, almost 10% of Bhutan’s glacierized area 
would vanish and the meltwater flux would drop by as much 
as 30%. Under the conservative scenario of an additional 
1°C regional warming, glacier retreat is going to continue 
until about 25% of Bhutan’s glacierized area will have dis- 
appeared and the annual meltwater flux, after an initial 
spike, would drop by as much as 65%. Citation: Rupper, S., 
J. M. Schaefer, L. K. Burgener, L. S. Koenig, K. Tsering, and E. R. 
Cook (2012), Sensitivity and response of Bhutanese glaciers to 
atmospheric warming, Geophys. Res. Lett., 39, L19503, 
doi: 10.1 029/20 12GL0530 10. 

1. Introduction 

[ 2 ] Glaciers are particularly sensitive to climate change, 
making them vulnerable elements of the environment. Of 
potential concern for societies is the rapid glacier retreat 
of Himalayan glaciers. Besides the vast potential for 
hydroelectric power, these glaciers contribute to the major 
rivers in Asia [Immerzeel et al., 2010; Kaser et al., 2006], 
and glacier meltwater is critical for agriculture in many, 
especially drier, regions in summer. Beyond, glacial lake 
outburst floods are one of the major natural hazards in this 
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region, and the hazard potential is rapidly growing as a result 
of glacier thinning and retreat [ Richardson and Raynolds, 
2000]. Nevertheless, estimates of ongoing and near- future 
glacier change in the Himalayas relative to the climate 
forcing remains poorly quantified and thus highly contro- 
versial [e.g., Cogley et al., 2010; Immerzeel et al., 2010; 
Jacob et al, 2012; Bolch et al., 2012]. 

[ 3 ] Recent studies have estimated large-scale Himalayan 
glacier change through remote sensing techniques. However, 
the results range from dramatic ice mass loss [e.g., Dyurgerov 
et al, 2009; Bolch et al., 201 1] with some estimates for high 
Asia as large as ~50 Gigatons per year [Matsuo and Heki, 
2010], to almost constant glacier mass [Jacob et al., 2012], 
These discrepancies could be partially explained by aspects 
inherent to differing remote sensing approaches, such as short 
observation periods, uncertainty in the remote sensing algo- 
rithms applied, and spatial resolution of the data. Detailed, 
regional glaciological mass-balance studies integrating multi- 
decadal glacier change are therefore desirable for comparison 
to both remote sensing studies as well as local, field-based, 
mass balance studies. One of the greatest challenges is the 
severe lack of field data for model and remote sensing vali- 
dation. This stems, in part, from the immense number of 
glaciers spread over a vast region, the complex politics of the 
region, and the rugged terrain. Modeling of glacier mass 
balance and sensitivity can complement, extend, and moti- 
vate field and remote sensing studies of mass balance. 
Therefore, quantifying the uncertainties that attend model 
estimates of glacier mass balance and sensitivity in the 
absence of accurate validation data is critical. From the per- 
spective of a glacier modeling approach to estimating glacier 
changes, these uncertainties stem from uncertainties in the 
glacier models, glacierized area, and climate data. Indeed, 
most studies use a single climate data set to quantify glacier 
changes and projections [e.g., Rupper and Roe, 2008; 
Immerzeel et al., 2010], without quantifying the uncertainties 
associated with these data. In this study, we assess the 
uncertainties in projected glacier changes associated with the 
climate data and derive a conservative scenario for recent and 
ongoing glacier change in the Kingdom of Bhutan, in the 
monsoonal Himalayas. 

[ 4 ] Bhutan is chosen for several reasons. First, Bhutan 
exemplifies an area where little data on glacier changes are 
available and where it is logistically difficult to obtain field- 
based studies, a common problem for many regions of the 
Himalayas. Few glaciological studies have addressed glacier 
changes in Bhutan, with the exceptions of a glacier area 
inventory [Mool et al., 2001], quantitative estimates of gla- 
cier retreat and area decrease from 1963 to 1993 for a sam- 
pling of-15% of Bhutan’s glaciers [Karma et al., 2003] and 
an atlas of glaciers of Bhutan [Iwata, 2010]. Not all of these 
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Figure 1. Equilibrium Line Altitude Sensitivity to Climate Change. All Panels: Black contours are elevations in 1000 m 
intervals. White outline is the Kingdom of Bhutan, (top) Average annual precipitation from the Climate Research Unit 
Climatology 2.0 (CRU CL 2.0) gridded data. Note the extremely high precipitation rates centered over Bhutan and adjacent 
regions, (middle) Change in glacier equilibrium line altitude (ELA) for a 1°C increase in temperature. Note the high glacier 
sensitivity to temperature change throughout the non-arid Himalayas (>~500 rnm/yr), peaking over Bhutan and the neighbor- 
ing Himalayas of Nepal, India, and Myanmar, (bottom) Change in glacier ELAs for an increase in precipitation equal to 


500 nun. Note that ELA sensitivity to even this large increase 

studies are peer-reviewed, and the uncertainties associated 
with the reports are not well known. All of these studies and 
reports do provide important insights into the rapid changes 
occurring in the large glacierized areas in Bhutan over the 
past half century. However, to date there are no published 
field- or modeling-based mass balance estimates for these 
glaciers, no estimates of glacier sensitivity to climate in the 
region, and no estimates of related change in meltwater flux 
over time. Thus, while there is significant evidence that the 
glaciers in Bhutan have been retreating over the past half 


in precipitation is quite low over the non-arid regions. 

century, we know little about the causes of those glacier 
changes and what the future glacier change will be. 

[5] Second, glaciers in Bhutan, just as neighboring gla- 
ciers in India, Nepal, and Southwest China, sit in the bulls- 
eye of high snow accumulation glaciers (Figure 1, top) 
(see Text SI for methods detail)/ Sensitivity tests using 
a temperature-melt model (Text SI) support prior work 


'Auxiliary materials are available in the HTML. doi:10.1029/ 
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that show that high accumulation regions are extremely 
temperature-sensitive [e.g., Rupper and Roe, 2008; Fujita 
and Nuimura, 2011]. Therefore, Bhutan’s glaciers form a 
highly suitable natural laboratory to investigate glacier sen- 
sitivity and response to temperature change in the mon- 
soonal Himalaya (Figure 1). 

[6] Finally, there are socio-economic reasons to focus 
work on the glacierized regions of Bhutan. The Kingdom of 
Bhutan is a country at the forefront of climate change miti- 
gation strategies (Bhutan was the first developing country to 
receive climate mitigation monies from the UN’s Least 
Developed Countries fund), and faces many of the economic 
and hazard challenges associated with glacier changes in 
the larger Himalayan region [Nayar, 2009]. For example, 
Bhutan’s main economic export is hydroelectric power and 
the powerplant safety and viability depends critically on 
prediction of glacial lake outburst floods and related miti- 
gation measures [Belding and Vokso, 2011], Therefore, the 
well-being of Bhutan’s society requires accurate estimates of 
glacier change and meltwater production. 

[ 7 ] This work provides a comprehensive study of glacier 
systems in Bhutan, and the first associated quantification of 
present and near-future glacierized area loss and related 
changes in melt water flux with particular emphasis on the 
sensitivity of these quantifications to the choice of different 
climate input data sets. We base our investigations on new 
glacier mapping from satellite imagery, a suite of gridded 
climate data, and a robust glacier melt model. This approach 
provides a range in mass balance estimates over the entire 
glacierized region. 

2. Results and Conclusions 

2.1. Glaciated Area 

[8] In the Himalayas, widely used glacier data sources 
are the World Glacier Inventory (WGI) [World Glacier 
Monitoring Service and National Snow and Ice Data Center, 
2011], the Global Land Ice Measurements from Space 
(GLIMS) [Armstrong et al., 2012], and the Natural Earth data 
(NE) [Patterson and Kelso, 2011], We merged these datasets 
into a single glacierized area, and then refined this merged 
area by using 1999-2010, June through September, Landsat 
Enhanced Thematic Mapper Plus (ETM+) images and 1999 
i-cubed eSAT imagery (see Text SI for additional details). 
The total, refined glacierized area for Bhutan and the Bhutan 
watershed is presented in Figure 2a, resulting in 1930 km 2 
and 3120 km 2 , ±3%, respectively. Our area estimates are 
larger than earlier estimates for Bhutan’s glacierized area 
[Mool et al., 2001], with the greatest discrepancy occurring 
in the perennial snowpack and glacier accumulation areas. 

2.2. Mass Balance and Glacierized Area Changes 

[ 9 ] In order to calculate the net mass balance for Bhutan’s 
glacierized area, we determine both, the melt and accumu- 
lation rates across the glacierized area. Here we use a com- 
monly applied degree-day melt model [Hock, 2003] to 
estimate melt rates across the glacierized region of Bhutan, 
modified to capture spatial variations in surface conditions. 
The model assumes melt is proportional to the sum of all 
temperatures greater than zero (the degree-days), and that the 
proportionality constant relating melt to the degree days 
varies depending on the glacier surface (e.g., dirty ice, clean 
ice, fresh snow). While the melt factors are uncertain without 


validation, the melt factors used here are likely conservative 
estimates. (See Text SI for a more detailed description of the 
melt model.) Total annual accumulation is calculated as the 
sum of all precipitation. The precipitation magnitudes and 
patterns in this region of complex topography are highly 
uncertain. Here we err on the side of maximizing accumula- 
tion by assuming all precipitation falls as snow. This mini- 
mizes the net imbalance across the region for the given 
precipitation data. 

[ 10 ] If the integrated mass balance is negative (positive), 
the glacierized area will decrease (increase); if glaciers are in 
steady state with climate the integrated mass balance equals 
zero. Mass balance for each 250 m x 250 m grid cell is 
calculated as the total incoming accumulation minus the 
total outgoing melt, averaged over the period 1980 to 2000 
(Figure 2c) (see Text SI and Table SI). This temporal range 
provides the average climate over which most Bhutanese 
glaciers have already adjusted to or are currently responding 
to, given estimates of glacier response times (see Text SI). 
In addition, defining the mean climatology using twenty 
years of data helps insure that a few anomalously warm 
years do not significantly inflate the results. 

[ 11 ] Errors in the mapped areas, modeling approach, and 
climate input data can all lead to considerable uncertainties in 
glacier meltwater flux and mass balance estimates. Here we 
present conservative values for glacier mass balance and, in 
turn, minimum values for glacier change in Bhutan by using 
the following estimates: First, our high elevation perennial 
snowpacks and glacier accumulation areas are larger than 
those of earlier glacierized area estimates [Mool et al, 2001], 
The use of Landsat imagery to map these regions likely 
increases the mapping accuracy as compared to the coarse 
resolution data used previously by Mool et al. [2001] (30 m 
grid resolution versus 1:50,000 scale topographic maps), 
but the differences in mapped regions may also reflect 
uncertainties in glacierized area. Importantly, these high and 
cold regions contribute very little to the total melt, but 
increase the net accumulation across the region significantly. 
Therefore, using upper bounds for the glacierized area esti- 
mate in higher accumulation regions minimizes the resulting 
glacier change. We assume the total glacierized area is rea- 
sonable for this study. Future work should focus on quanti- 
fying the uncertainties associated for these areas further. 
Second, we assume all precipitation falls as snow. This errs 
on the side of maximizing snow accumulation and mass 
balance. Third, we use conservative melt factors in relating 
temperature to melt in order to err on the side of minimizing 
melt and glacier change. Indeed, our modeled results for 
much of the region fall within the range of mass balance and 
melt rate estimates for similar glaciers in similar settings in 
the Himalayas [e.g., Bolch et al., 2011; Fujita and Nuimura, 
2011; W. Yang et al, 201 1]. Here we focus on the regional- 
averaged changes and, therefore, rely less on providing the 
accuracy necessary to confidently assess the mass balance of 
any single glacier. This approach to the mapping and mod- 
eling provides conservative mass balance estimates for the 
region of Bhutan. One exception is over glacierized regions 
with thick debris cover. Thin debris-cover lowers the albedo 
and increases melt rates, while the insolating effects of thick 
debris-cover cm) can reduce melt rates below that of clean ice 
[Scherler et al, 201 1]. The model used here only accounts 
for the influence of debris on albedo, and, therefore melt may 
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Figure 2. Mass balance of glacierized areas in the Bhutanese watershed, (a) Index map showing location of Bhutan relative 
to greater Himalaya, (b) Modeled average 1980-2000 annual mass balance (total annual snowfall minus total annual melt) 
using the CRU TS 3.1 gridded climate data as input. Note nonlinear mass balance color scale. Base map for Figures 2b and 
2c are a 7.5 arc-second digital elevation model (GMTED2010) (darker to lighter grays = lower to higher elevations). Bhutan 
is outlined by the yellow line (recently disputed borders dashed) and all watersheds that feed Bhutan are outlined by the red 
line, (c) Calculated changes in meltwater flux (red) and glacierized area (blue) as a result of step-wise temperature changes 
from 0 to 6°C. Bold lines are the ensemble model means (Tables SI and S2). The thin lines are the results using the CRU TS 
3.1 data, illustrated spatially in Figure 2c. The gray shading represents the spread in projected June, July, August warming 
between 1980 to 1999 and 2080 to 2099 (based on the MMD-A1B models) for South Asia [Christensen et al., 2007]. 


be overestimated in regions of thick debris cover (Text SI). 
We quantify the potential impact of this effect below. 

[ 12 ] Important additional uncertainty arises from errors in 
the climate input data, which is the focus of this study. In the 
absence of a dense network of weather stations, we use 
gridded climate data from Climate Research Unit Time Series 
3.1 (CRU TS 3.1) [Jones and Harris, 2011], WorldClim 
[Hijmans et al., 2005], and Aphrodite [Yatagai et al., 2009] 
to provide a suite of possible estimates of regional climate 
across the glacierized region (see Table S2 and Text SI). 
These are commonly used datasets in high Asia. The com- 
bination of possible climate estimates used here provides 
some test of the uncertainties associated with choice of cli- 
mate data used. 


[ 13 ] For the given mass balance model and assumptions 
outlined above, the area-averaged, net mass balance is 
— 1.4 ± 0.6 m/yr (Table SI). The lower range of the mass 
balance estimates is within the range of estimates from other 
mass balance studies in the Himalayas [e.g., Bolch et al., 
2012]. Importantly, this implies that the Bhutanese glacier- 
ized area is considerably out of balance with the recent cli- 
matology, or in other words, even if no additional warming 
were to occur, the glacierized area of Bhutan will decrease in 
order to reach steady state (mass balance equal to zero) with 
mean climate. We calculate the change in area required to 
reach steady state with recent mean climate and the current 
mean mass imbalance by stepwise removing the 250 m x 
250 m glacierized grids that have the most negative mass 
balance until the regionally integrated mass balance is equal 
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to zero (see Text SI for more details). Consistent with the 
approach above, this provides again the minimum change of 
Bhutan’s glacierized area required to reach steady state with 
modem mean climate, yielding a decrease of 40 ± 1 9% over 
the next decades, in the absence of any additional climate 
change (Figure 2c and Table SI). In steady state, this mini- 
mum area loss results in an annual glacial meltwater flux 
decrease by 71 ± 20%. The 95% confidence intervals pro- 
vided for the mean changes in mass balance, glacierized area, 
and melt volume flux (on the order of 25-50%) are due only 
to the differences in climate datasets, illustrating the sensi- 
tivity of these results to the uncertainties in climate data. 
However, even the most conservative values (Figure 2 and 
Table SI) show that glaciers in Bhutan are out of equilibrium 
and that changes in glacierized area and meltwater flux are 
going to be significant. 

[ 14 ] In order for the climate to compensate for even the 
most conservative scenario, precipitation would need to 
almost double, which is outside of most reconstructions from 
paleoclimate proxies or from future and past climate model 
estimates [e.g., Jiang et al., 2011 ;Hu etal., 2008; Hewitt and 
Mitchell, 1996; Johns et al., 2003; Christensen et al., 2007]. 
Alternatively, regional cooling would need to occur, which 
would contradict recent and current climate trends and 
atmospheric C0 2 increases [ Christensen etal., 2007; A Yang 
etal., 201 1], Therefore, the current glacier-climate imbalance 
will near-certainly result in a large reduction in glacierized 
area and melt water resources in Bhutan. 

[ 15 ] Given the fact that melt may be over-estimated for 
heavily debris covered areas (see Text SI), we repeat the 
mass balance, percent area change, and meltwater flux 
changes for only the clean glaciers for comparison. Clean 
glacierized areas are identified using the Landsat imagery 
and the model applied to only those areas. The most con- 
servative area-averaged net mass balance for only clean 
glaciers is —0.3 m/yr. In order to reach steady state with this 
mass imbalance, area must decrease by nearly 6% and 
annual melt water flux would decrease by at least 25%. The 
clean glacier scenario is likely a best case scenario, and still 
represents significant imbalance between the glacierized 
area and climate. Thus the results for the entire glacierized 
area and clean glaciers both indicate that, even if no addi- 
tional warming were to occur, the glacierized area of Bhutan 
must decrease significantly to reach steady state with recent 
mean climate. 

2.3. Future Climate and Glacier Scenarios 

[ 16 ] To test the sensitivity of the glacierized area to ongo- 
ing warming, we calculate the area and meltwater volume 
change for each step-wise 0.5°C increase in temperature up 
to 6°C (Figure 2 and Table SI) (see Text SI) for the full 
glacierized area. Even for a conservative warming of 1°C, 
at minimum 25% of the glacierized area will be lost and 
present meltwater flux will decrease, after an initial melt 
spike, to 35% of today’s value. A warming of 2.5°C, the 
average projected temperature change for the next century 
over South Asia [Christensen et al., 2007], results in the loss 
of more than half of the glacierized area. Importantly, annual 
meltwater flux becomes negligible. 

[ 17 ] As for the time-scale of change, most future climate 
scenarios project 1°C warming within the next 50 years 
[Christensen et al., 2007]. By this scenario, the estimates of 


glacier area and meltwater changes in Bhutan are likely going 
to occur within decades, putting considerable time pressure 
on acquiring accurate climate data and developing robust 
glacier models to improve the accuracy of the results pre- 
sented here. It also highlights the need to clearly identify 
regions where high efficiency projects to mitigate and adapt 
to glacier change in the monsoonal Himalayas are needed. 

3. Discussion 

[is] The lion’s share of the glacial melt in the monsoonal 
Himalayas occurs during the peak summer months, after the 
spring snowpack melt and before the peak torrential rains 
of the summer monsoon, so during a period where the 
regional water supply may depend predominantly on glacier 
melt. In some regions, the glacial melt creates excess flow 
in summer, contributing to annual river flooding hazards. 
In drier regions of the monsoonal Himalayas, decreases in 
meltwater may lead to significant drying of rivers during 
summer months [Immerzeel et al., 2010; Bolch et al., 2012], 
In consequence, any venture relying on water during peak 
summer months, including hydroelectric power generation 
or agriculture, will have to take the potentially large decrease 
in annual melt water flux into account. 

[ 19 ] The approaches and results presented here likely carry 
relevance for the wider monsoonal Himalayas where glacier 
sensitivity to climate change is similar (Figures 1, middle, 
and 1, bottom). 

[20] Despite the geographically and topographically com- 
plex distribution of glaciers throughout South Asia, the 
conservative scenarios of glacierized area and meltwater flux 
changes presented here are unlikely to be limited to the 
glaciers of Bhutan, because numerous glacierized regions of 
the monsoonal Himalaya are situated in similar geologic, 
glaciologic and climatic settings. In turn, our results strongly 
suggest that the conclusions of a recently published remote 
sensing study by Jacob et al. [2012], implying quasi-constant 
volume of large ice-masses (>100 km 2 ) in the Himalayas 
over the last decade, are highly unlikely to be representative 
of glacier changes in the monsoonal Himalayas in the near 
future. 

[21] Importantly, the large uncertainty that attends any 
given climate dataset represents a significant difficulty in 
quantifying glacier mass balance and associated glacier 
changes in the Himalayas. Even if the uncertainties in gla- 
cierized area and the mass balance model were minimized, 
uncertainty in the climate inputs would represent a significant 
uncertainty. This is true not only for gridded climate data, but 
also for weather station data that are extrapolated to glacier- 
ized areas. Continued focus on quantifying regional glacier 
mass balance and climate through the wider Himalaya region 
is called for to further reduce the uncertainties in glacier 
sensitivity to climate change. 
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